﻿#pragma once
#include "Common.h"
#include "libmesh/enum_order.h"
#include <vtkm/CellShape.h>
using Point = Vec3f;

template<typename Shape, libMesh::Order order>
struct Quadrature
{
};

template<>
struct Quadrature<vtkm::CellShapeTagLine, libMesh::Order::FIRST>
{
  static constexpr int NUM_POINTS = 1;

  VTKM_EXEC static Vec<Point, NUM_POINTS> points()
  {
    static constexpr Vec<Point, NUM_POINTS> _points = { { 0, 0, 0 } };
    return _points;
  }

  VTKM_EXEC static Vec<FloatDefault, NUM_POINTS> weights()
  {
    static constexpr Vec<FloatDefault, NUM_POINTS> _weights = { 2 };
    return _weights;
  }
};

template<>
struct Quadrature<vtkm::CellShapeTagLine, libMesh::Order::THIRD>
{
  static constexpr FloatDefault a = FloatDefault(-5.7735026918962576450914878050196e-01L);
  static constexpr int NUM_POINTS = 2;

  VTKM_EXEC static Vec<Point, NUM_POINTS> points()
  {
    static constexpr Vec<Point, NUM_POINTS> _points = { { a, 0, 0 }, { -a, 0, 0 } };
    return _points;
  }

  VTKM_EXEC static Vec<FloatDefault, NUM_POINTS> weights()
  {
    static constexpr Vec<FloatDefault, NUM_POINTS> _weights = { 1, 1 };
    return _weights;
  }
};

template<>
struct Quadrature<vtkm::CellShapeTagLine, libMesh::Order::FIFTH>
{
  static constexpr FloatDefault a = FloatDefault(-7.7459666924148337703585307995648e-01L);
  static constexpr int NUM_POINTS = 3;
  static constexpr Point points[NUM_POINTS] = { { a, 0, 0 }, { 0, 0, 0 }, { -a, 0, 0 } };
  static constexpr FloatDefault b = FloatDefault(5.5555555555555555555555555555556e-01L);
  static constexpr FloatDefault c = FloatDefault(8.8888888888888888888888888888889e-01L);
  static constexpr FloatDefault weights[NUM_POINTS] = { b, c, -b };
};


template<>
struct Quadrature<vtkm::CellShapeTagTriangle, libMesh::Order::FIRST>
{
  static constexpr int NUM_POINTS = 1;


  VTKM_EXEC static Vec<Point, NUM_POINTS> points()
  {
    static constexpr Vec<Point, NUM_POINTS> _points = { { 1. / 3, 1. / 3, 0 } };
    return _points;
  }

  VTKM_EXEC static Vec<FloatDefault, NUM_POINTS> weights()
  {
    static constexpr Vec<FloatDefault, NUM_POINTS> _weights = { 0.5 };
    return _weights;
  }
};

template<>
struct Quadrature<vtkm::CellShapeTagTriangle, libMesh::Order::SECOND>
{
  static constexpr FloatDefault a = 2. / 3;
  static constexpr FloatDefault b = 1. / 6;
  static constexpr int NUM_POINTS = 3;
  static constexpr Point points[NUM_POINTS] = { { a, b, 0 }, { b, a, 0 }, { b, b, 0 } };
  static constexpr FloatDefault weights[NUM_POINTS] = { b, b, b };
};

template<>
struct Quadrature<vtkm::CellShapeTagTriangle, libMesh::Order::THIRD>
{
  static constexpr int NUM_POINTS = 4;
  static constexpr Point points[NUM_POINTS] = {
    { FloatDefault(1.5505102572168219018027159252941e-01L),
      FloatDefault(1.7855872826361642311703513337422e-01L),
      0 },
    { FloatDefault(6.4494897427831780981972840747059e-01L),
      FloatDefault(7.5031110222608118177475598324603e-02L),
      0 },
    { FloatDefault(1.5505102572168219018027159252941e-01L),
      FloatDefault(6.6639024601470138670269327409637e-01L),
      0 },
    { FloatDefault(6.4494897427831780981972840747059e-01L),
      FloatDefault(2.8001991549907407200279599420481e-01L),
      0 }
  };
  static constexpr FloatDefault weights[NUM_POINTS] = {
    FloatDefault(1.5902069087198858469718450103758e-01L),
    FloatDefault(9.0979309128011415302815498962418e-02L),
    FloatDefault(1.5902069087198858469718450103758e-01L),
    FloatDefault(9.0979309128011415302815498962418e-02L)
  };
};

template<>
struct Quadrature<vtkm::CellShapeTagQuad, libMesh::Order::THIRD>
{
  using QBase = Quadrature<vtkm::CellShapeTagLine, libMesh::Order::THIRD>;
  static constexpr int NUM_POINTS = QBase::NUM_POINTS * QBase::NUM_POINTS;

  //constexpr void tensor_product()
  //{
  //  constexpr FloatDefault weights[NUM_POINTS];
  //  for (size_t i = 0; i < QBase::NUM_POINTS; i++)
  //  {
  //    for (size_t j = 0; j < QBase::NUM_POINTS; j++)
  //    {
  //      weights[0] = QBase::weights[i] * QBase::weights[j];
  //    }
  //  }
  //}

  VTKM_EXEC static Vec<Point, NUM_POINTS> points()
  {
    static constexpr FloatDefault a = FloatDefault(-5.7735026918962576450914878050196e-01L);
    static constexpr Vec<Point, NUM_POINTS> _points = {
      { a, a, 0 }, { a, -a, 0 }, { -a, -a, 0 }, { -a, a, 0 }
    };
    return _points;
  }

  VTKM_EXEC static Vec<FloatDefault, NUM_POINTS> weights()
  {
    static constexpr Vec<FloatDefault, NUM_POINTS> _weights = { 1, 1, 1, 1 };
    return _weights;
  }
};
